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Abstract. In our previous work [1], we have proposed two methods for computing the 
luminosity distance d/^ in ACDM model. In this paper, two effective quadrature algorithms, 
known as Romberg Integration and composite Gaussian Quadrature, are presented to cal- 
culate the luminosity distance d^^^ in the Chevallier-Polarski-Linder parametrization(CPL) 
model. By comparing both the efficiency and accuracy of the two algorithms, we find that 
the second is more promising. Moreover, we develop another strategy adapted for approxi- 
mating in flat ACDM universe. To some extent, our methods can make contributions to 
the recent numerical stimulation for the investigation of dark energy cosmology. 
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1 Introduction 

The computation and numerical evaluation of distances is frequently encountered in the 
research of cosmological phenomena. In practice, it is common to compute the various 
cosmological distances as a function of the redshift z under certain cosmological models. 
Current cosmological observations indicate that the expansion of universe is accelerating, 
has a prominent dark energy content ACDM ~ 0.7 and is spatially flat [2]. Further, a 
Lambda cold dark matter (ACDM) model fits the data well, and is frequently used as a 
fiducial or background model. In the ACDM model, various cosmological distances can be 
expressed in terms of the elliptic integrals [3, 4]. 

Some works have focused on the luminosity distance d_L in the ACDM model and de- 
rived numerical approximations for the efficient and accurate evaluation of dL{z) given the 
cosmological parameters Q^n and $7a [5-7]. The computation of is useful in the analysis of 
distance-redshift relations of type la supernovae, and the approximation for di can also be 
directly used in the evaluation of other distances, for instance the angular diameter distance 
dA_ or the comoving distance r [4]. 

Recent years, Chevallier, Polarski [8] and hinder [10] proposed a simple parametrization 
of the dark energy equation of state (known as CPL): 



which is involved in the luminosity distance. The best fit values of wq and Wa are -1.58 and 
3.30[11], respectively. CPL parametrization is widely applied into both observational and 
theoretical analysis, and it has the talent to test the dynamics of many dark energy models. 
More discussions about CPL model can be seen [9-13]. In this paper, we just concentrate 
on the numerical analysis of the complicated integral contained in the luminosity distance 
of CPL parametrization model. Note that analogous integral may be encountered in many 
cases, e.g. the dynamical age of the universe or the angular diameter distance. 

For the rest of this paper we will pay our main attention to perform the numerical 
quadrature algorithms. Section 2 is a brief review of the luminosity distance in the CPL 
parametrization model. In section 3, we present two different quadrature algorithms and 



■w{z) =1^0 + 11), 



z 



(1.1) 
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compare their efficiency and accuracy based on the personal computer. Another approximate 
recipe of the luminosity distance in ACDM is developed in section 4. Finally, we discuss some 
improvements of the algorithms and possible extensions to other cosmological models briefly. 



2 Luminosity Distance in CPL Model 

In order to study the different dark energy models, the widely used method is assume an 
ad hoc equation of state w{z) = px/px for dark energy and parametrize w{z) [14]. CPL 
parametrization model was first proposed by M. Chevallier, D. Polarski [8] and E. V. Linder 
[10], and the parameterized w{z) can be written as equation 1.1. Thus the dark energy 
density px is given by 

pxiz) = p\f{z) (2.1) 

with 

/(,) = (1 + ,)3{i+«.o+^.) exp(-^). (2.2) 

For a spatially flat universe (A; = 0), the Friedmann equation can be expressed as 

H\z) = ^-^{pM + Px) = HlE\z) 

= Hl[^^{l + zf + VL^{z)l (2.3) 

where H{z) is the Hubble parameter, Vl^ is the dark matter parameter, VLj^{z) represents the 
time-dependent dark energy parameter, E[z) is the expansion rate of the universe. 

With the continuous equation and Friedmann equation 2.3, ^\{z) can be deduced as 

^k{z) = UAfiz), (2.4) 

where Q\ is the dark energy parameter at present time and f{z) is defined as equation 2.2. 
Hence, the luminosity distance in the CPL model can be written in the form 

,CPL _ c(i + z) r dz' 

The expression of the luminosity distance is complicate and there is no critically analytical 
solution for general parameter choice {wQ,Wa)- Meanwhile, we can take the luminosity dis- 
tance in ACDM universe as a special case with [wq^Wo) = (—1.0,0.0). From this point of 
view, the will degenerate to 

. c{l + z) r dz' 



dl = ' ' / . (2.6) 

^0 Jo V^m(l + z'f + 

3 Quadrature Algorithms 

Equation 2.5 is just a one-dimensional integral, but there are four variables in the integrand, 
i.e. z, u)o, Wa- Because of the different variations of the three variables {yi^,WQ,Wa)-, 
approximating the integral directly seems to be impossible if we want to obtain desirable 
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accuracy. On the other hand, if we approximate the integrand with multivariables inter- 
polation technique [15] and then integrate the approximate polynomial, the expression will 
remain complicate and is hardly to implement in practice. 

Instead of developing an analytical approximation, an effective quadrature algorithm 
may be more helpful. We present two conventional numerical integration methods, known as 
Romberg Integration and Gaussian Quadrature [16] , in the following subsections and compare 
their performances to see which one is more suitable for calculating the luminosity distance. 

For simplifying the following description, we just consider the integral in equation 2.5, 
defined as : 



fE 



dz' 

W) 



(3.1) 



with 



,3(l+'W0+«'a) 



exp(- 



'l + z 



where fim + Vt\ = 1. 



3.1 Romberg Integration 

Romberg Integration gives preliminary approximations with the Composite Trapezoidal rule 
and then applies the Richardson extrapolation process to improve the accuracy. For each 
integer k = 2,3,4, ...,n and j = 2,3, ...,/c, an 0{h^j^) approximation formula can be written 
as 



f{x)dx = Rk,, + 0{h?^), 



where /i^ = (6 — a) ^ and the iterative formula R^j is 



■k,j 



R 



+ 



4''-i - 1 ■ 



(3.2) 



(3.3) 



In order to obtain the complete iterative process, we should considerate the case when j = 1- 
In general, the R^^i is provided by using the trapezoidal approximation, then we have 



Ri,i = ^[f{a) + f{b)], 



and 



2k- 



Rk-1,1 + hk-i f{a + {2i - l)hk) 



i=l 



(3.4) 



(3.5) 



for A; = 2, 3, n. 

The main results generated from the above formulas are listed in Table 1. The Romberg 
Integration has an additional desirable feature that it allows an entire new row in the table 
to be calculated by performing one additional application of the Composite Trapezoidal rule. 
Then it uses an averaging of the previously calculated values to obtain the succeeding entries 
in the row. The method used to construct a table of this type calculates the entries row by 
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Table 1. Approximation results for Romberg Integration. 



^1,1 










-^2,1 


R2,2 








-^3,1 


R3,2 


-^3,3 






-^4,1 


Ra,2 


-^4,3 


Ri,A 




Rn,l 


Rn,2 


Rn,3 


Rn,A 


Rn,n 



row, that is, in the order R2,i, R2,2, etc. R. L. Burden describes a detailed algorithm 
in [17]. 

We can preset an integer n to determine the number of rows. In many cases, however, it 
is confused to ensure whether the output is satisfactory or too many entries are unnecessary 
to generate. For using the iterative technique sufficiently and saving the running time, 
we can set an error tolerance for the approximation and generate n, within some upper 
bound, until some consecutive entries agree to within the tolerance. In this paper, we choose 
\Rn,n — Rn,n~i\ < to generate the approximations. 

The Romberg Integration method can be used in conjunction with other numerical 
quadrature formulae to obtain successive improved values. Further, the method is applicable 
to a very large class of functions. 



3.2 Composite Gaussian Quadrature 

Gaussian Quadrature chooses the points for evaluation in an optimal, rather than equally 
spaced. The nodes ri,r2,...,rn in the interval [a,b] and coefficients ci, C2, c^, are chosen 
to minimize the expectancy obtained in the approximation 

rb " 

/ fix)dx^^cJir,). (3.6) 
i=i 

With the Legendre polynomials, we can determine the nodes and coefficients easily. The 
nodes are the zeros of the nth Legendre polynomial Table 2 lists the nodes and 

coefficients for n=7 and 8. More detailed calculations and values of nodes and coefficients 
can be found in [16, 18]. Such quadrature formulae is called Gauss-Legendre formulae. 

The error term of higher-order derivative in quadrature formulae of higher degree is 
difficult to evaluated or even boundless, the quadrature formulae 3.6 is not recommended 
to obtain desirable accuracy, although the Gaussian Quadrature is stable. Instead, we can 
divide the interval [a,b] into some subintervals [xi,Xi+i], and apply the low-order Gaussian 
Quadrature to each subinterval. Then, summing over all the values as the final approxima- 
tion, we have 

/ fix)dx = E / ^(^)'^^' (3-^) 

where a = xq < 2:1 < ... < x^. = b and the subscript m donates the number of the 
subintervals. 
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Table 2. Nodes and weight coefBcicnts for Gauss-Legendre integration. 



n 


nodes r.„^j 


Coefficients Cn,j 


7 ± 


0.9491079123 


0.1294849662 


± 


0.7415311856 


0.2797053915 


± 


0.4058451514 


0.3818300505 







0.4179591837 


8 ± 


0.9602898565 


0.1012285363 




0.7966664774 


0.2223810345 


± 


0.5255324099 


0.3137066459 




0.1834346425 


0.3626837834 



Let the subintervals be of equal size, the composite Gaussian Quadrature can be written 



as: 



f{x)dx = 2 Z]' 



m— 1 



2 2 ^ 



(3.8) 



where h = {b — a)/m donates the size of the subinterval. 

The Gaussian Quadrature formulae can be applied only when /(x) is explicitly known, 
so that f{x) can be evaluated for any desired value of x. Naturally, orthogonal polynomial 
other than the Legendre polynomials also can be used, such as the Gauss - Chebyshev, Gauss 
- Jacobi and Gauss - Hermite formulae. 



3.3 Performance of the Tv^o Algorithms 

In the section, we perform the efficiency and accuracy of the two quadrature algorithms with 
Fortran program. We create a sample containing about 10^ redshift data which ranges from 
to 1100 to present a quantitative comparison. 

The Romberg Integration is based on equation 3.3, and we choose the error tolerance 
\Rn,n — Rn,n-i\ < 10~^ to generate the approximation. With the fast convergence rate of the 
iterative, there is no obvious distinctness, including the execution time, if we set the error 
tolerance \Rn,n — Rn,n-~i\ < 10"^ instead. Different parameters {m,wo,Wa) are chosen to 
evaluate the composite Gaussian Quadrature formulae 3.8, but we fix the number of nodes 
n = 8. 

Fig 1 shows the relative error of different parameter choices for Gaussian Quadrature 
which illustrates that more subinterval division can improve the accuracy of the algorithm and 
the relative error also depends on the choice of wq and Wa- However, the error is insensitive 
to the m for redshift 2: < 50. Actually, the general Gaussian Quadrature based on equation 
3.6, i.e. m = 1, is precise enough to evaluate the integral values in this case. Therefore, in 
order to obtain desirable accuracy and efficiency, we can set a greater number n and suitable 
m to extend the composite Gaussian Quadrature to a wider redshift distribution. 

The main results of the two algorithms are listed in Table 3. Note that the Gaussian 
Quadrature obviously takes less time than Romberg Integration if the same accuracy is 
required. However, extra interpretation about Table 3 should be emphasized. The execution 
time and efficiency just reflects the relative results of the codes with each other, and depends 
on the compiler as well as the different computing environment. More discussions are available 
in [1] and [6]. 
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Figure 1. The absolute relative error as a function of z for different parameter choices (m^wi^^Wa)- 
The error is no more than 0.5% when to > 80 for the best fit parameters u>o — —1.58, Wa ~ 3.3, = 
0.3. The black dashed line, i.e. (P, -1.0, 0.0), denotes the accuracy of the polynomial approximation 
(see section 4). 



Table 3. The main results of the Romberg Integration, Gaussian Quadrature and polynomial ap- 
proximation. The last two columns show that the iterative of Romberg Integration spends more time, 
although its accuracy is well under control. 



type 


parameters 


time(s) 


maximum error (%) 


Romberg Integration 


(-1.58, 3.3) 


210.960 


0.04 


Gaussian Quadrature 


(80, -1.58, 3.3) 


7.457 


0.48 




(100, -1.58, 3.3) 


9.266 


0.23 




(150, -1.58, 3.3) 


13.822 


0.028 




(200, -1.58, 3.3) 


18.564 


0.003 




(30, -1.00, 0.0) 


1.888 


0.43 




(50, -1.00, 0.0) 


3.151 


0.16 


Polynomial approximation 


(P, -1.0, 0.0) 


0.016 


0.41 



4 A Polynomial Approximation to 

The general analytical expression of the luminosity distance in the spatially flat ACDM 
model is given by equation 2.6. Because it is frequently used in practice, many papers have 
focused on the numerical analysis of the integral equation. In this section, another polynomial 
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approximation is described for the similar considerations. Following the notation introduced 
in our previous work [1], substituting s = y^(l — fl^) /^Ira and u = 1/z' into equation 2.6 
yields 



4 



l + z 



where 



c/Hq y/sW^ 



T{t) 



du 



^/u"^ + u 



(4.1) 



However, we note that the behavior of T(r) has some deficiencies. First, the derivative 
of T(r) becomes singular as r — t- 0"*". Second, the domain of T{t) extends to infinity. Either 
one is detrimental to the approximation using polynomials. Fortunately, such shortages can 
be eliminated by proper change of variables. Based on the consideration above, we introduce 
a mathematical transformation a; = l/(r + l) further to constraint the variation of x within 
< X < 1, which leads to 



l + z 



c/Hq a/sH^ 



/( 



1 



where 



JO 



s/{l + z) + l' 



du 



fi 



1 



s + l 



(4.2) 



\/(l - n) • (3^2 - 3n + 1) 

The main purpose of the section is to obtain the approximate expression of the equation. 

Since the integrand intends to infinity as n — )• 1~, we utilize a roundabout strategy to 
achieve the suitable polynomial. By analyzing the integrand, we find that it can be factorized 
as above (equation 4.2), which inspires us to write it as: 



fix) 



El 



(4.3) 



where Ojit* is a polynomial defined as: ajU* 

One prominent aspect of the equation 4.3 is that it can be integrated analytically. For 
quarrying out the six best-fitting free parameters Oj, we must impose some constraints upon 
the equation 4.3. With the modish range of Q^n in mind, we just minimize the relative error 
between function / and / with 0.1 < < 1- The restrictions what we adopt are listed as 
following: 



fil) 



fi 



1, 



'3' 



/(I) 
/(I) 



/(I), 



(4.4) 



where 1/3 is approximately equal to the minimum value of x within the considered parameter 
space, i.e., x G [1/3,1]. 

Utilizing equation 4.4, we can obtain a linear equations which contains just two free 
parameters. For instance, we can choice 04 and 05 as the variables to be determined. In 
order to derive the total parameters, we define the relative error as: 



_ 4_ 



(4.5) 
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Figure 2. Approximate and analytical function / {upper panel) and the absolute percentage relative 
error {lower panel). The solid and dashed lines represent the / and / with 0.1 < < 1 in the upper 
panel, respectively. 



Just as [5] has pointed out, the error tends to be dominated by z — )■ 0, which yields 

e(x)maa; = 1^ - lU^O- (4.6) 

The square of e(x), namely e^(x), as a continuous function, is more convenient for us to 
acquire the relation between the rest two parameters 04 and 05. From the mathematical 
theorem we know that the first derivative of e^(x), which contains three variables x, and 
05, must be equal to zero strictly if it reaches the local maximum. By solving the three 
nonlinear equations, the most appropriate relation between 04 and 05 can be written as 
following: 

04 = 12.15722 - 2.92471 • 05, (4.7) 

where as is still unknown to us. Because we are more interested in the minimum of the 
e{x)max-: substituting equation 4.7 into equation 4.6 and regulating the valve of 05 make us 
find the most desirable outcome of e{x)„iax 

to be about 0.96% when ag = -44.63290. 
However, we have emphasized that the implicit requirement of the coincidence of func- 
tion values at end points could be unnecessarily strong [1]. If we control the relative error to 
the minimum when ilm is taken as the best observational value 0.3, the final approximation 
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log z 



Figure 3. The distribution of the global relative error in with different z and ilm- The maximum 
error, dominated by the small redshift (z < 0.1), is less than 1%. 



we construct can be expressed as: 

f{x) = y/l-x ■ (S.llSOTx^ - 22.69338x^ + 21.09474x3 

-6.03039x2 _^ 0.32109X - 2.80713) + 2.80713, (4.8) 

where 

1 3 1 ~ 

s/{i + z) + r ^ 

Fig. 2 shows the results of the fit for the function / (see equation 4.2) and its residual, which 
is not more than 0.15%. 

With the definition of the relative error in equation 4.5, the distribution of the global 
relative error based on the polynomial approximation (equation 4.8) with various redshift and 
Jim is plotted in Fig. 3. As we can see from the figure, the maximum error, which is less than 
1%, tends to be contributed by the small redshift z. Similarly, we compute the running time 
and accuracy of the approximate luminosity distance based on the same criterions presented 
in section 3.3. From Table 3 we conclude that the polynomial spends less time than the other 
two methods, though its error can't be improved further. However, Fig. 1 shows that the 
accuracy of the polynomial at low redshift is inferior to the composite Gaussian Quadrature. 
And as the increase of redshift z, the relative error tends to be around 0.09%. 
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Because approximating the integrand directly can decrease the error of the target- 
integral, one promising extension of our method is to the linear growth factor 5{a) [19, 20] 
contained only fim and The general form of 6{a) in ACDM universe can be written as 



where a and H{a) represent the scale factor and Hubble parameter, respectively. Kasai [7] has 
developed an effective recipe for evaluating it recently. However, a more compact form may 
be helpful to research the matter density perturbation with numerical simulation efficiently, 
and as a candidate our idea behind the method may be useful. 

5 Conclusions 

Two different quadrature algorithms are presented to evaluate the integral involved in lumi- 
nosity distance in CPL parametrization model. Through the comparison of the efficiency and 
accuracy, the composite Gaussian Quadrature is more promising to apply into the evaluation 
of such complex integral. Because of the generalization of the algorithms, we can also extend 
them to other parametric models of dark energy, for instance the two-index parameterizations 
by Huterer et al.[22], and the four- index parameterizations by Hannestad et al.[23]. 

Additionally, because general Gaussian Quadrature is precise enough for refshift z < 50, 
it is unnecessary to divide the interval [a, b] to be of equal size. The adaptive integration 
combined with the composite Gaussian Quadrature may be more helpful. 
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